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Abstract 

We compare theoretical and simulation results for static and dynamic properties for a model 
of particles interacting via a spherically symmetric repulsive ramp potential. The model displays 
anomalies similar to those found in liquid water, namely, expansion upon cooling and an increase 
of diffusivity upon compression. In particular, we calculate the phase diagram from the simulation 
and successfully compare it with the phase diagram obtained using the Rogers- Young (RY) closure 
for the Ornstein-Zernike equation. Both the theoretical and the numerical calculations confirm the 
presence of a line of isobaric density maxima, and lines of compressibility minima and maxima. 
Indirect evidence of a liquid-liquid critical point is found. Dynamic properties also show anoma- 
lies. Along constant temperature paths, as the density increases, the dynamics alternates several 
times between slowing down and speeding up, and we associate this behavior with the progres- 
sive structuring and de-structuring of the liquid. Finally we confirm that mode coupling theory 
successfully predicts the non-monotonic behavior of dynamics and the presence of multiple glass 
phases, providing strong evidence that structure (the only input of mode coupling theory) controls 
dynamics. 
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I. INTRODUCTION 



Water and some other liquids exhibit anomalous behavior close to their freezing lines 
Their phase diagrams have regions characterized by a negative thermal expansion 
coefficient, i.e., these liquids expand upon cooling at certain temperatures and pressures. 
Besides the density anomaly, such liquids also have other peculiar thermodynamic and dy- 
namic behaviors |3[ . For example, the isothermal compressibility increases upon cooling and 
the diffusivity increases upon pressurizing. Usually the region of the diffusion anomaly is 
wider than the region of the density anomaly, so that the latter is completely contained in 
the former ^|. The anomalous behavior of the thermodynamic properties of water has been 
connected to the existence of a hypothetical liquid-liquid critical point in deeply supercooled 
states 0, 0, In the case of water, this critical point is located in an experimentally 
unaccessible region. Recently, using the potential energy landscape formalism, it has been 
argued [5] that under certain assumptions on the statistical properties of the potential energy 
landscape, the existence of a density anomaly must lead to the existence of a liquid-liquid 
critical point. 

In the case of water, anomalies are thought to be related to the tetrahedrality of the in- 
terparticle potential. On average, each water molecule has four nearest neighbors connected 
by hydrogen bonds. However, tetrahedrality is not a necessary condition for anomalous be- 
havior and several spherically symmetric potentials that are indeed able to generate density 

and/or diffusion anomalies have been proposed 0, H H H 0, Q H Q • 

Interestingly, 

such potentials may be purely repulsive, providing evidence that different microscopic mech- 
anisms can generate density anomalies. These potentials can be regarded as the simplest 
models which yield water-type thermodynamic and dynamic anomalies and it is important 
to fully characterize their thermodynamic and dynamic behavior. An additional advantage 
in studying spherical potentials is that their behavior can also be studied within a theoreti- 
cal framework, for their thermodynamic properties can be calculated using accurate integral 
equation closures. 

Here we study, using extensive molecular dynamics (MD) simulations, a specific, spheri- 
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12. 13. 141 with the aim of 



cally symmetric repulsive potential introduced by Jagla 
fully characterizing both static and dynamic extreme loci in the temperature-density plane. 
We complement the MD study with integral theory calculations based on the (thermody- 
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namically consistent) Rogers- Young (RY) closure, which is known to give accurate results 



for other repu 



potential 
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sive potentials, such as the square shoulder potential 



171 ] and the star polymer 



19l l20j . We also compare numerical results in the low T region with predic 
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which we solve using the RY 



tions based on the ideal mode-coupling theory (MCT) 
static structure factors as input. We find that MCT is able to predict the non-monotonic 
behavior of dynamics and the presence of multiple glass phases, providing further evidence 
that structure (the only input of MCT) controls dynamics in this system. 

II. METHODS 



A. Discrete Molecular Dynamic Simulations 

We study the linear ramp potential introduced by Jagla \lQ 
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00 r < ctq 

U(r)= { U Q {a x - r)/(<7i - O ) ^0 < r < 01, (1) 
r > 0i 

focusing on the specific choice of A = 01/00 = 1.76, which has been studied previously by 
Jagla |lOj . Using Monte Carlo simulations, Jagla showed that a density maximum is found 
in the liquid phase. For this choice of A, the potential in Eq. (JTJ is a good candidate for 
studying the connection between thermodynamic and dynamic quantities. In order to study 
the dynamic properties, we apply the discrete MD method, approximating the continuous 
potential by a sequence of step functions with n small vertical steps, 



UJr) 



00 r < O 

kAU a 1 -(k + |)Ar <r <a x -{k- \)Ar 
r > 0i - 



2> 

At 
2 



(2) 



where Ar = (01 - O )/(n + 1/2), AU = U /(n + 1/2), and k = 1,2,3,. . .n. The unit of 
length is 00, while Uq is the unit of energy. Temperature is measured in units of energy, i.e., 
k B = 1. Simulation time is measured in units of cr ^/ m /Uo, with m as the particle mass, 
and pressure in units of Uo/<Jq. The number density is defined as p = N/L 3 , where L is the 
size of the simulation box and N is the number of particles. 
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The standard discrete molecular dynamics algorithm has been implemented for particles 



interacting with step potentials [23, 



25, 26]. ong straight 

lines with constant velocities. When the distance between the particles becomes equal to 
the distance for which U(r) has a discontinuity, the velocities of the interacting particles 
instantaneously change. The algorithm calculates the shortest collision time in the system 
and propagates the trajectories of particles from one collision to the next. Calculations of 
the next collision time are optimized by dividing the system into small sub-systems, so that 
collision times are computed only between particles in neighboring sub-systems. Since the 
total energy E is rigorously conserved, it is best to study the NVE-ensemble in the cubic 
box of a fixed volume V = L 3 with periodic boundary conditions. 

We consider N = 1728 particles in our simulation. For constant temperature and pressure 
simulations, the Berendsen thermostat and barostat are used. Configurations in the NVE 
ensemble are saved on a logarithmic spacing for further processing, after discarding an initial 
equilibration period for a time larger than the correlation time at the particular state point 
[27| . The diffusion coefficient, measured in units of aQ^/Uo/m, is calculated as 

(K , + Q-rW V 

where (. . denotes an average over all particles and over all t'. The dynamic structure 
factor for a given vector q is defined as 



5(q,t) = /^^e^(iq-[r i (t' + t)-r i (0])\ , (4) 

where (. . .) t > denotes the average over all t' and S(q) = S(q, 0) is the static structure factor. 
The normalized structure factor 

is called the density correlator. The isotropy of the liquid allows us to average 5*(q, t) over 
different q with the same modulus. In the following, we bin together all q within a mesh 
Aq = ir/L. 



B. Rogers- Young Closure 



Integral equation theories are powerful tools for studying the structure and thermody- 



namic properties of liquids 



28l |29|. One assumes a two-body interaction potential for the 
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particles and introduces the total pair correlation function h(r), related to the pair distribu- 
tion function g(r) = h(r) + 1, and the direct correlation function c(r). The goal is to solve 
the self-consistent Ornstein-Zernike equation, 

h{r) = c(r) + p J dr'c{\r — r'\)h(r'), (6) 

where both c(r) and h(r) — or, equivalently, g(r) — are unknown. Equation (jUJ) in Fourier 
space takes the form 

h(q) = c(q) + pc(q)h(q) } (7) 

where c(q) is related to S(q) as c(q) = (1 — l/S(q))/p. 

According to the particular form of the interaction potential U(r), one can choose a 
certain ansatz for c(r), which relates it to the interaction potential and to h(r) and allows 
one to solve Eq. (jUJ) analytically in some cases j3] and numerically otherwise. The two 
frequently employed closures, Percus-Yevick (PY) or the hypernetted chain (HNC) [29], 
suffer from being thermodynamically inconsistent 28j]. This means that although they 
can provide good estimates of the static structure factor, they cannot be reliably used for 
determining the phase diagram of the system. More sophisticated closures have thus been 
developed in recent years, which have a built-in thermodynamic consistency. This is achieved 
by introducing an extra parameter in the ansatz, which can then be determined to fulfill 
such a condition. 

The RY closure jsjj belongs to the thermodynamically consistent group of closures ob- 
tained by appropriately mixing PY and HNC through the parameter (. The ansatz for c(r) 
in terms of h(r) becomes 



c(r) = exp[—(3U(r)] 



1 exp[(h( r) - c(r))/(r)] - 1 



where f(r) is the 'mixing function', 



[h(r)-c(r) + l], (8) 



/(r) = l-exp(-Cr). (9) 

For C = 0, one recovers PY closure, while for ( —>■ oo Eq. (jHJ reduces to the HNC condition. 

The OZ equation can be solved using Eq. (jHJ) for a given value of (. The correct solution 
corresponds to that value of ( for which the compressibilities Kt, calculated using the 
"virial" and the "fluctuations" routes agree, ensuring thermodynamic consistency. This 
allows us to reliably use the RY closure not only to calculate the static structure factor, but 



also the phase diagram, as we will do in the following. Moreover, RY gives particularly good 
results for a purely repulsive potential, such as the studied ramp potential. It has already 
been successfully tested against simulations for square shoulder potential s Il7 ], and for both 
experiments and simulations for the star polymer effective potential BET 



C. Mode Coupling Theory 

The density correlator defined in Eq. (jHJ) is the fundamental quantity of interest in MCT, 
a set of generalized coupled Langevin equations which can be closed within certain approx- 
imations |2ll . I22I ] . Interesting behavior of these observables arise when the dynamics of the 
system becomes slower, i.e., when the dynamical behavior is of the "supercooled" type j^. 
A typical two-step relaxation occurs for the density correlators on approaching dynamical 
arrest, indicating the emergence of two distinct time scales in the system's structural re- 
laxation. A first relaxation process, the f3 relaxation, occurs at short times, and is due to 
particles exploring the cages formed by their nearest neighbors. A second relaxation, the a 
relaxation, occurs at longer time scales, when particles are able to escape the cages. MCT 
predicts the existence of a glass transition at a characteristic temperature Tmct, where 
the time scale of this second relaxation r a diverges so that the particles will always remain 
trapped in their cages. For T < Tmct the correlators do not relax any more, reaching a finite 
plateau value at long t, defined as the nonergodicity parameter f(q) = lim^oo 0(g, t); f(q) 
jumps discontinuously from zero to a finite (critical) value f c {q) at the transition, signaling 
the occurrence of an ergodic (fluid) to a nonergodic (glass) transition. The transition is 



kinetic, i.e., nothing happens at the thermodynamic properties of the system c 
MCT predictions are often found to be in agreement with experimental 



ose to T MC T- 
I and simula- 



tion results [34], although in real systems the a-relaxation time does not diverge, but only 
becomes increasingly larger. This is due to the intervention of other processes, commonly 
termed "hopping" processes, which restore ergodicity and are not included in the MCT 
treatment of the ideal glass transition described above. 

In mathematical terms, the non-ergodicity parameters f(q) are the long time solutions 
of the MCT equations, i.e., 

/((/) -m(q), (10) 



1 " /(?) 
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where the memory kernel m(q) is quadratic in the correlator 

1 C rl 3 k 

m(q) = -J ^pV(q, k)/(*)/(|q - k|), (11) 

where k — |k|. The vertex functions V are the coupling constants of the theory, which are 
given only in terms of the static structure factor and number density of the system, 

V(q, k) = L [ q . ( q - k ) c(|q - k|) + q • k c(k)} 2 x 

S(q)S(k)S(\q-k\). (12) 

Equations (fTU|) and (fTT|) define a system of non linear equations with a trivial solution 
f(q) = 0. However, for certain values of the vertex functions, the solutions have a bifurcation 
point, locating the glass transition. At this transition point a solution f(q) > emerges. 
The time evolution of the density correlators is found by solving the full MCT equations, 

4>(q, t) + fi 2 (g)0(g, t) + [ m(q, t - t')<j>(q, t') dt' = , (13) 

Jo 

where Q 2 (q) = q 2 /((3mS(q)) and m(q,t) is the time-dependent memory kernel 

1 f d 3 k 

m(q,t) = - J — V( q ,k)<KM)<K|q-k|,t). (14) 

The two-step relaxation is well described by MCT through an asymptotic study of the 
correlators near the ideal glass solutions. The a-relaxation is effectively described by a 
stretched exponential, 

<P(q,t) = A q exp[-(t/r q )% (15) 

where the amplitude A q determines the plateau value and f3 q < 1. 

MCT predicts a power-law divergence of the a-relaxation time as well as a power-law 
decrease of the diffusivity as the system approaches the ideal glass transition. In this study, 
we compare the dynamic behavior evaluated from the MD simulations with corresponding 
MCT predictions. 

III. RESULTS 

A. Dependence on the Number of Steps in the Ramp 

We first study how the results of molecular dynamic simulations for the potential U n {r) 
in Eq. (J2J) converge to the results for the continuous ramp potential U{r) in (JTJ), as n — > oo 
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where n is the number of approximating steps. The exact equation of state for systems 
described by a pair potential can be calculated from the partition function, i.e., by the 
integral f f ... f exp[— $^i>j Ufa — r^/ksT] YliLi d> r i- Replacing the continuous potential 
U(ri — Tj) by a step-function is analogous to replacing an integral by a sum in a rectangular 
approximation which is known to converge to the integral as 1/n 2 . Thus we can expect that 
the pressure P n obtained in the discrete molecular dynamics for given n converges to the 
value for the continuous potential as 1/n 2 . On the other hand, the probability for a 
particle to jump over a step of size AU is proportional to exp(—AU/k B T) ~ 1 — Uo/nksT. 
The diffusion coefficient must be a different iable function of this probability. Thus we can 
expect that the diffusion constant D n in discrete MD approaches its limiting value as 
1/n. Figures Ufa) and|2^b) confirm these predictions. In the following, we limit ourselves to 
the case n = 144 which, as shown in Figs. Efa) andEfb), is sufficiently close to the n —>■ oo 
case. 

B. Structure Factors and Comparison with RY 

Figure H3 shows the density dependence of the structure factor, as calculated from the MD 
at T = 0.063. At this low T the liquid, even at low densities, is significantly structured, as 
shown by the large amplitude of the first peak. On increasing the density, contrary to the 
normal liquid behavior, the amplitude of the first peak is reduced. As discussed below, the 
destructuring of the liquid is associated to a speed-up of the particle dynamics. If p is further 
increased, the second peak becomes dominant and (as discussed in the following) its increase 
correlates with a slowing down of the dynamics. The first peak is significantly reduced and 
acts as a pre-peak on the major peak. The crossover of the leading amplitude from the first 
to the second peak resembles the behavior previously discovered in star-polymers of large 
functionality |3J| . In these systems, the star-star interaction can be effectively reduced to an 
ultra-soft repulsive logarithmic interaction |l8j]. However, the peak positions of S(q) change 
positions with density in this opposed to what is found in ramp potential. 

Figures E[a)-BJd) compare the static structure factor calculated using simulation data 
along with the results of the RY-closure for several densities at T = 0.063. A tolerance of 
5.10 -5 is used for the thermodynamic consistency. The RY integral theory provides a correct 
description of the p dependence of the dynamics. RY correctly predicts the structure factor 
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at high densities, while for intermediate and low densities RY gives a fairly good estimate of 
the structure factor except that the first peak in the structure factor is lower than the one 
calculated from simulation. RY thus tends to underestimate the structure of the liquid. 



C. Phase Diagram 

Figure El shows the phase diagram of the model obtained by slow cooling (thin lines 
representing isochores) and by accurate simulations of individual state points (circles). The 
simulations for each density and temperature are initialized with random configurations and 
are equilibrated to the desired temperature using Berendsen's thermostat for a sufficient 
period of time. The equilibration time was estimated as the time when the density correlator 
0(q, t) at the first peak of S(q) decays to zero. After equilibration, each configuration was 
left to run at constant energy for a time dependent on the speed of the dynamics, covering 
at least ten equilibration times. The temperature of maximal density (TMD), shown by 
a bold line, was obtained by connecting the points on each isochore corresponding to the 
minimal pressure, since these points correspond to the points of minimal volume due to the 
general relation 



and (dV/dP)T < for a mechanically stable system. Thus, volume increases upon cooling 
at constant pressure in the region to the left of the TMD line. 

At low T, the system spontaneously crystallizes during the time of the simulation in 
different crystalline forms. Consistent with Jagla's calculation (Q|), we find that: 

(i) At low densities, the system crystallizes into a face centered cubic crystal structure 
upon slow cooling, when the temperature reaches the values indicated by crosses. The 
crystallization is marked by a sharp drop in the potential energy, associated with a 
fast release of latent heat. The line connecting these points can be regarded as a 
line of homogeneous nucleation. It has a marked negative slope, corresponding to the 
smaller density of the crystalline phase compared to liquid. The equilibrium melting 
line, which also has negative slope, is located at much higher temperatures, so that 
a large portion of the density anomaly lies in the supercooled state. This situation 
is completely analogous to the situation in water, in which the region of the density 




(16) 
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anomaly is located near the negatively-sloped freezing line. In the ramp model, the 
anomaly does not exist below p = 0.2432, where the line of the minimal density merges 
with the line of the maximal density. Interestingly, the slope of the homogeneous 
nucleation line becomes positive below this density, as in a normal liquid where the 
crystal phase has a larger density than in the liquid phase. This behavior may exist 
in stretched water at negative pressures but has never been observed experimentally 
or in simulations. 

(ii) In the region of intermediate densities 0.317 < p < 0.352 the liquid does not crystallize 
and enters the glassy state below T < 0.042. We are not able to equilibrate the system 
below this temperature. 

(hi) For p > 0.352, the system crystallizes into a rhombohedral crystal after being equili- 
brated for a long time at T = 0.041. Interestingly, the density anomaly vanishes at 
densities slightly lower than the density at which this new crystal phase emerges. This 
crystalline phase is characterized by parallel columns formed by equidistantly spaced 
atoms. The spacing among these atoms is slightly larger than the hard core distance. 
Thus in this crystalline phase each atom has two neighbors in its repulsive ramp. The 
projection of these columns onto a perpendicular plane forms a triangular lattice with 
spacing approximately equal to the diameter of the repulsive ramp. The columns are 
shifted with respect to each other by one third of the nearest-neighbor spacing so that 
the atoms form three crystalline planes, perpendicular to the columns. In each of these 
planes atoms form a triangular lattice with a spacing a/3 times larger than the spacing 
in the triangular lattice formed by the projection of the columns. 

(iv) A third distinct (hexagonal) crystalline phase is observed for p m 0.7865 (outside 
the range of densities explored in Fig. 0) . At this state point the liquid, after some 
initial equilibration, crystallizes into a crystalline phase characterized by a hexagonal 
symmetry of one of its crystalline planes. This crystalline type is also formed by 
parallel columns consisting of equidistantly spaced atoms. The spacing among these 
atoms is slightly larger than the hard core distance. The projection of these columns 
onto a perpendicular plane forms a hexagonal lattice with spacing slightly smaller than 
the hard core. This spacing is equal to v3/2 of the distance between atoms in the 
columns. The neighboring columns are shifted with respect to each other by one half 
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of the atom spacing in the columns, so that the atom and its two neighbors in the 
neighboring column form an equilateral triangle. Thus in this crystalline type, each 
atom has eight nearest neighbors in its repulsive ramp: two in the same column and 
six in the three neighboring columns. 

(v) For even larger densities not studied in this work, hexagonal close packed and finally 
hard-sphere face cubic centered crystals are expected [10]. 

We now discuss the possibility of a liquid-liquid critical point in this model. By quadrat- 
ically extrapolating the isochores into the glassy state, we observe a crossing at very low 
temperatures. This crossing of the near density isochores is equivalent to dP/dV\T = 
and hence predicts the existence of a critical point with coordinates T pa 0.025, P ~ 0.838, 
and p pa 0.346, close to the largest density at which an isochore develops a negative slope. 
Interestingly, the TMD line, if extrapolated, appears to go directly to this putative critical 
point. 

Another characteristic feature is the behavior of the isothermal compressibility, which 
shows an anomalous increase upon cooling between the lines of maximal and minimal com- 
pressibility, are shown in red in Fig.^Ja). The part of the low density branch with a positive 
slope corresponds to a compressibility minima, while the high density branch corresponds to 
the compressibility maxima. As in water and other materials, the isothermal compressibility 
line crosses the TMD line at the point of its maximal temperature (vertical slope) due to the 
mathematical properties of the second derivative of the equation of state j^fj. Again, as in 
water, the line of compressibility maxima, if extrapolated, seems to approach the putative 
critical point. 

Since the RY closure is thermodynamically consistent, it is possible to evaluate the phase 
diagram theoretically as well. Indeed, once the static structure factor S(k) is obtained, the 
pair correlation function g(r) can be calculated by taking the inverse Fourier transform of 
S(k). All other thermodynamic quantities are calculated using the inter-particle potential 
function U(r) and g(r). In Fig. 03(b) we show the phase-diagram consisting of eight isochores 
of the system calculated using RY. Comparing Fig. [31(a) and Fig. 03(b), we can see that the 
RY predicts a smaller value of the pressure and also that the density maxima points are 
shifted to higher volumes. Despite this discrepancy in P, the shape of the isochores is very 
well reproduced. Even in the RY case, density maxima line appears in the phase diagram. 
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While in the simulation, extremely slow dynamics prevents access to the region where the 
critical point is probably located, in the RY calculations, no convergence of the parameter 
( is achieved in the same region. Again, a smooth extrapolation of the calculated isochores 
is consistent with a crossing point, and hence a critical point. 

The thermodynamic behavior discussed above is analogous to the behavior of the one- 
dimensional model for which an exact solution can be obtained (see Appendix). 



D. Dynamics 

Next we focus on the dynamic properties of the model. Figure Efa) shows the lines of 
diffusivity maxima and minima, i.e., the locus of points where dD/dP\r = in the phase 
diagram. The region of the anomalous increase of D upon compression (dD/dP\T > 0) 
embraces the region of density anomaly as it does in two-dimensional models as well as in 
water [JQ. 

Figure |Hf a) shows the density dependence of D along several isotherms in a density range 
extending to p = 0.7. While only one minimum and one maximum are observed at high T, 
at low T, data suggest the possibility of a more complex behavior of the density dependence 
of D. Indeed, for the range of densities 0.492 < p < 0.579, the diffusion coefficient drops 
sharply for T < 0.067, and then recovers at higher densities, suggesting the presence of more 
than one locus of diffusion maxima. It is interesting to observe that this second maximum 
separates regions with different underlying crystal structures (orthorhombic and hexagonal) 
(37 1 , similar to the location of the low p maximum, separating regions where the fee and 
orthorhombic crystals are respectively stable. 

It is interesting to compare the p dependence of the characteristic times numerically 
calculated with the prediction of mode coupling theory using RY results for S(q) (MCT- 
RY). A comparison with the dynamics predicted by MCT can be performed by comparing 
the T and p dependence of D with the T and p dependence of the (inverse) characteristic 
time scale of decay of density fluctuations. This time can be calculated within MCT, from 
the solutions of Eq. (fT3"|) . More specifically, for each wavevector q, the decay time can be 
defined as the value at which <f)(q,t) reaches the value 1/e. Figure E^b) shows the inverse 
relaxation time at the wavevector corresponding to the first peak of the structure factor r-f 1 
as a function of p for different T. The same sequence of minima and maxima characteristic 
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of the MD diffusivity data is also found in the prediction of MCT-RY. This agreement is 
consistent with the possibility that the structure factor, the only input in the MCT, also 
encodes the system's dynamic properties. 

To better characterize the low temperature dynamics and investigate the possibility of 
different glasses, we next study the decay of the density autocorrelation functions <p(q, t) 
as a function of t at T = 0.063 in Fig. [7J and calculate <f){q,t) from the MD trajectories 
Fig. Ufa). The non-monotonic behavior of D and T\ discussed in the previous figure is also 
seen in the decay of 4>(q,t). At the highest studied density p > 0.787, the correlator does 
not decay to zero in the time scale of our simulations. For lower 0, the decay becomes faster 
until = 0.352 and then it starts to slow down again. For < 0.26, crystallization prevents 
the approach to a glass state. The decay of the correlation functions (at q corresponding 
to the first peak of the structure factor) can be well represented via a stretched exponential 
decay [Eq. (fi3j) ] at very high densities, p = 0.702 (/3 = 0.89), and by simple exponential 
decay {(3 = 1) at low p. A similar behavior is seen in Fig. Efb) where the predictions of 
the MCT equations are shown. Again, on decreasing p, first a speed-up and then a slowing 
down of the dynamics is observed. 

To quantify the comparison between MD and MCT-RY, we show the density dependence 
of T\ and T2 in Fig. |H] (the characteristic time of the first and second peaks of the structure 
factor respectively) both from MD and MCT-RY. The two sets of data show the same 
anomalous behavior, showing the minimum corresponding to the maximal diffusivity at the 
same p. 

It is interesting to observe that while the slowing down of dynamics can be numerically 
followed for a large dynamical range at a high density, a low density crystallization prevents 
the generation of a glass structure. Unlike in MD, RY solutions can be found in a wider 
density range. It is thus interesting to ask whether a glass line is predicted by MCT-RY 
at low densities and, if so, how the two glasses differ. By solving the MCT equations for a 
wider range of densities, we find two distinct glasses at T = 0.052: p = 0.682 and p = 0.257. 
These two glasses are characterized by different critical non-ergodicity parameters (shown in 
Fig. EJ). While the low density glass non ergodicity factor is similar to the one found in hard 
sphere systems, the high density glass is characterized by large amplitudes both at the first 
and second peak of S(q). This resembles the one found in star polymers at high density, 
though the amplitude of the second peak in the latter becomes larger than that of the first 
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peak. 



IV. SUMMARY AND CONCLUSIONS 

In this work we have presented a systematic numerical study of the static and dynamic 
properties of a system of particles interacting with a spherical repulsive potential, with a 
range of interaction of the order of the particle diameter. The simplicity of the model 
makes it possible to study it with efficient numerical algorithms and, theoretically, with 
self-consistent integral theories |2.8i]. Our results accomplish the following: 

(i) They confirm the previous theoretical calculations by Jagla concerning the ex- 
istence of a line of density maxima in the phase diagram of this potential. Results 
also provide an accurate evaluation of this line as well as of the line of compressibility 
maxima and minima. 

(ii) They confirm that different crystal structures are found at low temperatures, depend- 
ing on the density Q], and provide estimates of the homogeneous nucleation line for 
the liquid-crystal transitions. 

(iii) They provide evidence of the possibility of a liquid-liquid critical point at T « 0.025, 
P ~ 0.838 and p ~ 0.346. The location of the critical point is below the homogeneous 
nucleation line or glass transition line and cannot be accessed by simulations. A the- 
oretical RY calculation of the phase diagram is able to reproduce the thermodynamic 
anomalies. These calculations also suggest the possibility of a liquid-liquid critical 
point but, again, its precise location cannot be determined due to the impossibility 
of equating the "virial" and "fluctuation" compressibilities with enough accuracy in 
this region of the phase diagram. Thus the existence of the critical point proposed by 
the extrapolation of the equation of state into the deeply supercooled region remains 
questionable j^. 

(iv) They show that the RY closure agrees reasonably well with simulations of the sys- 
tem with the repulsive ramp potential and that it reproduces the static and dynamic 
anomalies. The RY closure slightly underestimates the pressure, and shifts the anoma- 
lies to a region of lower density. 
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(v) They provide evidence that dynamic anomalies in this model have a structural origin 
and they are indeed captured by the MCT-RY theory, which uses only structural in- 
formation as input. Diffusion anomalies are encountered before the density anomalies, 
consistent with the case of water j^J. 

(vi) They suggest the possibility of different types of glasses in this simple system, consis- 
tent with the existence of different crystalline phases. A more extensive study based on 
the MCT-RY may help evaluate the location and the intersections (39, 40, 41 1 between 
different glass transition lines. 
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APPENDIX A: ONE DIMENSIONAL SYSTEM 



Applying the Takahashi method |42j, one can find the equation of state for the one- 
dimensional system of particles interacting via the ramp potential [see Eq. ([TJ]. In this case, 
the partition function \1/ for the Gibbs potential G(P,T) can be factored, so that 

G = ~NM(P,T), (Al) 

where 

POO 

#(P, T) = / exp (-p (U(r) + Pr)) dr. (A2) 
Jo 

Substituting Eq. ((TJ) into Eq. (jA2|) and integrating, we find that 

,w p r r i \ exp[-(3(Pa + Up)} - exp(-/3Pa 1 ) exp(-/? PaQ 

(3(P - Pc ) + JP ' {A6) 

where (3 = l/k B T and P c = C/ /(<ri - cr ). Since V = (dG/dP) T 

(dV/dP) 



P = N/V = - P* , (A4) 



which, after differentiation, leads to 



e -f3(U + a P) ^ _ 



F P -l3(Uo+v P)\—£o I I 1 1 P -p<nP\zx 21 i_ 1 I 1 v ; 

e HP-Pc) T /3(P-P C ) 2 J ^ e (P-Pc) T PP 2 /3(P-P C ) 2 \ 
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Note that P c plays the role of a critical pressure, at which the equation of state becomes 
discontinuous for T — > 0. Indeed, 



Hence, an infinitesimal increase of pressure around P c and small T leads to a finite increase 
of density from for P < P c to l/<7o for P > P c . Thus, in the one-dimensional system, 
there is an analog of a critical point at T = T c = 0, P = P c , and p = p c , at which the 
isothermal compressibility diverges. At P = P c , Eq. (jA5|) simplifies to 



For P < P c there exists a region in which there is a density anomaly and a region in which 
there is a compressibility anomaly, exactly as in the three-dimensional model studied in this 
work. 

The isochores described by Eq. (IA5J) are plotted in Fig. El One can see that the tem- 
perature of maximal density approaches the critical point as P — > P c . 
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FIG. 1: Schematic representation of the repulsive ramp potential Q and its discontinuous version 
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(a) 




1/n 

FIG. 2: (a) Pressure P as a function of n~ 2 [where n is denned in Eq. (|2"|)] for T = 0.063 and 
p = 0.260. (b) Diffusion coefficient D as a function of n _1 for the same state point. 
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FIG. 3: Density dependence of the structure factor at T = 0.063. Note the progressive reduction 
of the amplitude of the first peak and the progressive increase of the second peak on increasing p. 
Different curves refer to p = 0.272, 0.296, 0.322, 0.352, 0.384, 0.421, 0.464, 0.512, 0.567, 0.629, and 
0.702. 
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FIG. 4: Comparison of the static structure factor, obtained in simulations and in the RY closure 
for (a) low density p = 0.260, (b)-(c) intermediate densities p = 0.340 and p = 0.502, and (d) high 
density p = 0.702. 



23 



(a) MD 



0.8 



Q-0.6 



0.4 





TMD 




DM 




K T M 


X 


FCC 


■ 


RH 





— TMD 


p 


= 0.296 


<M P 


= 0.287 


P 


= 0.269 


>->P 


= 0.260 




= 0.252 


x — x p 


= 0.244 


*— * P 


= 0.237 


G-OP 


= 0.229 



FIG. 5: (a) Phase diagram of a system of particles interacting via the potential defined in Eq. (J2J) 
with n = 144. Lines indicate P{T) for several isochores at the following values of p, from top 
to bottom: 0.378, 0.364, 0.352, 0.340, 0.328, 0.317, 0.306, 0.296, 0.287, 0.277, 0.269, 0.260, 0.252, 
0.244, 0.242, 0.240, 0.237, 0.229, 0.223. Also shown are the locus of density maxima (8V/dT\ P = 0) 
(bold blue line), the locus of compressibility maxima and minima dV/dP\x = (dashed red 
line), and the locus dD/dP\x = (bold dotted line). At low p, lines terminate when the system 
crystallizes or does not equilibrate within the available simulation time. Note that different crystals, 
fee and rhombohedral, are found for p < 0.317 and p > 0.352, respectively. The extrapolated 
isochores at large p cross at a finite T (open box), consistent with the possibility of a liquid- 
liquid critical point for T below the ones that can be investigated numerically, (b) Corresponding 
phase diagram obtained within the RY closure. An extrapolation of the isochores also shows the 
possibility of a very low temperature liquid-liquid phase transition. 

24 



Q -2 
o 



O 



-6 



(a) 



MD 




*— * T=0.042 
<y^> T=0.066 
>^< T=0.063 
C3-Q T=0.067 
H — h T=0.068 
T=0.074 
T=0.084 
T=0.105 
T=0.126 
T=0.147 
T=0.168 



2 0.3 0.4 



MCT-RY 





0.052 


x— x T = 


0.063 


<-< T = 


0.073 


£^T = 


0.084 


^^T = 


0.105 


□— □ T = 


0.157 


G-GT = 


0.210 




0.8 



0.8 



FIG. 6: (a) The behavior of the diffusion coefficient as a function of density p for several isotherms 
obtained from MD simulations, (b) The density autocorrelation time Tf 1 from MCT-RY. At low 
temperatures when the density is very high or low, t\ increases sharply, showing the transition to 
a glassy state. Note the similar behavior of D and r-f 1 . 
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FIG. 7: The behavior of the correlators for the wave vector q\ corresponding to the first peak of 
the static structure factor for T = 0.063 for (a) MD and (b) MCT-RY. The dashed horizontal line 
indicates the 1/e value used to define the characteristic time t\. 
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FIG. 8: The density behavior of the correlation times (a) ri, and (b) T2, obtained directly from 
MD simulations for T = 0.063 and for the MCT calculations based on RY. The behavior of the 
diffusion coefficient is also shown. All curves show the anomalous decrease with density for small 
densities. 
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FIG. 9: The behavior for T = 0.052 of the critical non-ergodicity parameter f q for the two critical 
densities p = 0.257 and p = 0.682. 
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FIG. 10: The P — T phase diagram of the one-dimensional system with the ramp potential (ctq = 
l,o"i = 1.76, Uq = 1) defined in Eq. (Q). The thin solid lines are isochores for p = 0.05,0.1, ...,0.95 
from bottom to top. The bold dotted lines are isochores for p = \jo\ and p = 2/ (do + 0"i) 
between which the density anomaly region bounded by the temperature of the maximal density 
line (bold blue line) is located. The dashed bold line bounds the region of anomalous isothermal 
compressibility extrema (both minima and maxima). 
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